/*
 * =====================================================================================
 *
 *       Filename:  sa.cpp
 *
 *    Description:  Implementation of simulated annealing algorithm
 *
 *        Version:  1.0
 *        Created:  06/07/12 13:32:01
 *       Revision:  none
 *       Compiler:  gcc
 *
 *         Author:  Moos Hueting (mh), mooshueting@gmail.com
 *   Organization:  UCL
 *
 * =====================================================================================
 */
#include "sa.h"

#include <math.h>
#include <time.h>
#include <iostream>

#include "features/char_features.h"
#include "bivgauss/bivgauss.h"
#include "constraintSet.h"

#define CHAR_SIZE 17
#define START_TEMP 2000
#define TEMP_FACTOR 0.99

using namespace std;
using namespace cv;

Mat chomp(Mat src) {
    Size s = src.size();

    for (int i = 0; i < s.height; i++) {
        for (int j = 0; j < s.width; j++) {
            if (src.at<double>(i, j) < 0)
                src.at<double>(i, j) = 0;
            else if (src.at<double>(i, j) > 1)
                src.at<double>(i, j) = 1;
        }
    }

    return src;
}

SA::SA() {
    SA("hiragana.csv");
}

SA::SA(const char *csv) {
    // initial temperature
    temperature = START_TEMP;

    // seed rng
    srand(time(0));

    // allocate space for solution
    currentSolution = Mat(CHAR_SIZE, CHAR_SIZE, CV_64FC1, Scalar(0));
    currentIso = isocontour(currentSolution);

    // setup multivariate normal
    CvMLData ml;
    ml.read_csv(csv);
    Mat mData = ml.get_values();
    mData.convertTo(mData, CV_64F);
    MN = MultiNorm(mData.t());

# if 0
    RNG rng(rand());
    // randomly fill solution
    rng.fill(currentSolution, RNG::UNIFORM, 0, 0.6);
# endif

# if 0
    Mat bg = bivGauss(CHAR_SIZE, 0, 0, 3, 3, 0, 1.0);
    currentSolution = currentSolution + bg;
    bg = bivGauss(CHAR_SIZE, 0, 0, 1, 1, 0, 1.0);
    currentSolution = currentSolution - bg;
# endif

}

SA::SA(Mat _currentSolution, double _temp) {
    currentSolution = _currentSolution;
    currentIso = isocontour(currentSolution);
    temperature = _temp;

    // seed rng
    srand(time(0));
}

void SA::nsteps(int nSteps, int nTemp) {
    // perform n number of steps, decreasing
    // the temperature every tempStep steps

    for (int i = 0; i < nSteps; i++) {
        step();  // do step
        if (!(i % nTemp))
            lowerTemperature();  // lower temperature every nTemp steps
    }
}

void SA::step() {
    // generate new solution
    Mat newSolution = mutate().clone();
    Mat newIso = isocontour(newSolution).clone();

# if 0
    double diffEnergy = energy(newSolution) - energy();

    if (diffEnergy / energy() < 0.2)
        currentSolution = newSolution;
# else
    // find its energy
    double newEnergy = energy(newIso);

    double rValue = static_cast<double>(rand())/RAND_MAX;
    bool t = (exp((newEnergy - energy()) / temperature) > rValue);
//    cout << t << endl;

    // if the new energy is lower, keep regardless
    if (newEnergy >= energy()) {
        currentSolution = newSolution.clone();
        currentIso = newIso.clone();
    } else if (t) {
        // otherwise, use exponential function to
        // check if we want to keep new state anyway
        currentSolution = newSolution.clone();
        currentIso = newIso.clone();
    } else {
        currentIso = isocontour(currentSolution);  // WHY IS THIS NECESSARY?
    }
#endif
}

double SA::energy(Mat res) {
    // if character is empty, return default value
    if (countNonZero(res) == 0)
        return -1e12;

    // constraints work with 8UC1, so prepare
    Mat tmp = res * 255;
    tmp.convertTo(tmp, CV_8UC1);

    // low energy means good solution
    // return constraints.getError(tmp);
    int rCorners = corners(tmp);
    double rHeight = height(tmp);
    double rWidth = width(tmp);
    double rConvexity = convexity(tmp);
    double rComplexity = complexity(tmp);
    int rDiffSegs0 = diffsegs(tmp, 0);
    int rDiffSegs1 = diffsegs(tmp, 1);
//    int rAvgThickness0 = avgThickness(tmp, 0);
//    int rAvgThickness1 = avgThickness(tmp, 1);
    int rEndpoints = endpoints(tmp);
    double rRatio = ratio(tmp);
    int rHoles = holes(tmp);
    int rDC = diagConnections(tmp);
    int rCC = connectedComponents(tmp);
//    double rAxisMean0 = axismean(tmp, 0);
//    double rAxisMean1 = axismean(tmp, 1);
    double rDensity = density(tmp);
//    int rJunctions = junctions(tmp);
    double rVariance0 = variance(tmp, 0);
    double rVariance1 = variance(tmp, 1);
//    int rHoleNoise = holenoise(tmp);
    double rGradOrient0 = gradOrientation(tmp, 0);
    double rGradOrient1 = gradOrientation(tmp, 1);
    double rGradOrient2 = gradOrientation(tmp, 2);
    double rGradOrient3 = gradOrientation(tmp, 3);
    double rGradOrient4 = gradOrientation(tmp, 4);
    double rGradOrient5 = gradOrientation(tmp, 5);
    double rGradOrient6 = gradOrientation(tmp, 6);
    double rGradOrient7 = gradOrientation(tmp, 7);
    double rCurvature = curvature(tmp);

    // WARNING: ORDER MATTERS (must match order of CSV file)
    Mat x = (Mat_<double>(24, 1) << rCorners , rHeight, rWidth, rConvexity ,
                rComplexity , rDiffSegs0 , rDiffSegs1 ,
                rEndpoints , rRatio , rHoles ,
                rDC , rCC , rDensity , rVariance0 ,
                rVariance1 , rGradOrient0 , rGradOrient1 , rGradOrient2 ,
                rGradOrient3 , rGradOrient4 , rGradOrient5 , rGradOrient6 ,
                rGradOrient7 , rCurvature);

    // cout << x << endl;

    return MN.like(x);
}

Mat SA::mutate() {
      // add or subtract gauss
      // generate random mean value for gaussian (x and y)
      double mux = -(CHAR_SIZE - 1)/2 + ((CHAR_SIZE-1)/2) *
                        (static_cast<double>(rand())/RAND_MAX);
      double muy = -(CHAR_SIZE - 1)/2 + ((CHAR_SIZE-1)/2) *
                        (static_cast<double>(rand())/RAND_MAX);

      // generate random standard deviation
      double sx = ((CHAR_SIZE-1)/2)*(static_cast<double>(rand())/RAND_MAX);
      double sy = ((CHAR_SIZE-1)/2)*(static_cast<double>(rand())/RAND_MAX);

      // generate random correlation (between -1 and 1)
      double r = -1. + 2*static_cast<double>(rand())/RAND_MAX;

      // maximum value of gaussian
      double weight = 0.1 * static_cast<double>(rand())/RAND_MAX;

      Mat bg = bivGauss(CHAR_SIZE, mux, muy, sx, sy, r, weight);

      // decide if we subtract or add
      if (rand() % 2) {
          Mat newSol = currentSolution.clone() + bg;
          return chomp(newSol);
      } else {
          Mat newSol = currentSolution.clone() - bg;
          return chomp(newSol);
      }
}

void SA::lowerTemperature() {
# if 1
    temperature *= TEMP_FACTOR;
# else
    temperature -= 10000;
# endif
}
